
library(stargazer)
library(tidyverse)
library(lmtest) # for coeftest
library(sandwich) # for calculation of the cluster-robust se
library(patchwork)
library(Cairo)
library(stringr)
library(effects)
library(nnet) # for multinom
my_red  <- '#E83015' # 猩々緋
my_blue     <- '#005CAF' # 瑠璃
my_green    <- '#86C166' # 苗
my_yellow   <- '#F7C242' # 花葉
my_purple 	<- '#6A4C9C' # 桔梗

my_red  <- '#CB1B45'
my_blue     <- '#0089A7'

mlogit.clust <- function(model,data,variable) {
  beta <- c(t(coef(model)))
  vcov <- vcov(model)
  k <- length(beta)
  n <- nrow(data)
  max_lev <- length(model$lev)
  xmat <- model.matrix(model)
  # u is deviance residuals times model.matrix
  u <- lapply(2:max_lev, function(x)
    residuals(model, type = "response")[, x] * xmat)
  u <- do.call(cbind, u)
  m <- dim(table(data[,variable]))
  u.clust <- matrix(NA, nrow = m, ncol = k)
  fc <- factor(data[,variable])
  for (i in 1:k) {
    u.clust[, i] <- tapply(u[, i], fc, sum)
  }
  cl.vcov <- vcov %*% ((m / (m - 1)) * t(u.clust) %*% (u.clust)) %*% vcov
  return(cl.vcov = cl.vcov)
}

se_return   <- function(model, data){
    se_mat  <- summary(model)$standard.errors
    se_cluster  <- sqrt(diag(mlogit.clust(model, data, 'region_id')))
    se_mat[1, ]     <- se_cluster[1:(length(se_cluster) / 2)]
    se_mat[2, ]     <- se_cluster[(length(se_cluster) / 2 + 1):length(se_cluster)]
    return(se_mat)
}


suv_df  <- read.csv('governor_app_analysis_df_0521.csv') 
suv_df$gov_change_new3  <- relevel(factor(suv_df$gov_change_new2), ref = 'stay')
suv_df$gov_change_rem3  <- relevel(factor(suv_df$gov_change_rem2), ref = 'stay')
suv_df  <- suv_df %>% mutate(grp_logged = log(grp), democracy_lower = media_freedom	+ economic_liberalization + civil_society, elite_consolidation = political_pluralism + elites + political_institutions, factor_year = as.factor(year))


df_fit1     <- suv_df %>% select(gov_change_new3, democracy_lower,  growth, grp_logged, unemp_delta, term_limits_dum, nat_res_2006, population_2018, distance_from_moscow_km, rus_prop, age, tenure, outsider, gov_opp, region_id, factor_year, out_prop) %>% na.omit()
fit1_multinom   <- multinom(gov_change_new3 ~ democracy_lower + growth + grp_logged + unemp_delta + term_limits_dum + log(nat_res_2006 + 1) + log(distance_from_moscow_km + 1) +
              rus_prop + age + tenure + outsider + gov_opp + out_prop + factor_year, data = df_fit1)
df_fit2     <- suv_df %>% select(gov_change_new3, democracy_lower,  growth, grp_logged, unemp_delta, term_limits_dum, ur_share, ur_share_2003, nat_res_2006, population_2018, distance_from_moscow_km, rus_prop, age, tenure, outsider, gov_opp, region_id, factor_year, out_prop) %>% na.omit()
fit2_multinom   <- multinom(gov_change_new3 ~ democracy_lower + ur_share + ur_share_2003 + growth + grp_logged + unemp_delta + term_limits_dum +  + log(nat_res_2006 + 1) + log(distance_from_moscow_km + 1) +
              rus_prop + age + tenure + outsider + gov_opp + out_prop + factor_year, data = df_fit2)
df_fit3     <- suv_df %>% select(gov_change_new3, democracy_lower,  growth, grp_logged, unemp_delta, term_limits_dum, ur_fed, ur_share_2003, nat_res_2006, population_2018, distance_from_moscow_km, rus_prop, age, tenure, outsider, gov_opp, region_id, factor_year, out_prop) %>% na.omit()
fit3_multinom   <- multinom(gov_change_new3 ~ democracy_lower + ur_fed + ur_share_2003  + growth + grp_logged + unemp_delta + term_limits_dum +  + log(nat_res_2006 + 1) + log(distance_from_moscow_km + 1) +
              rus_prop + age + tenure + outsider + gov_opp + out_prop + factor_year, data = df_fit3)
df_fit4     <- suv_df %>% select(gov_change_new3,  growth, grp_logged, unemp_delta, term_limits_dum, last_margin, nat_res_2006, population_2018, distance_from_moscow_km, rus_prop, age, tenure, outsider, region_id, factor_year, out_prop) %>% na.omit()
fit4_multinom   <- multinom(gov_change_new3 ~  last_margin + growth + grp_logged + unemp_delta + term_limits_dum + log(nat_res_2006 + 1) + log(distance_from_moscow_km + 1) +
              rus_prop + age + tenure + outsider + out_prop , data = df_fit4)
df_fit5     <- suv_df %>% select(gov_change_new3,  growth, grp_logged, unemp_delta, term_limits_dum, margin_all, nat_res_2006, population_2018, distance_from_moscow_km, rus_prop, age, tenure, outsider, gov_opp, region_id, factor_year, out_prop) %>% na.omit()
fit5_multinom   <- multinom(gov_change_new3 ~ margin_all + growth + grp_logged + unemp_delta + term_limits_dum + log(nat_res_2006 + 1) + log(distance_from_moscow_km + 1) +
              rus_prop + age + tenure + outsider + gov_opp + out_prop + factor_year, data = df_fit5)
df_fit6     <- suv_df %>% select(gov_change_new3, elite_consolidation,  growth, grp_logged, unemp_delta, term_limits_dum, nat_res_2006, population_2018, distance_from_moscow_km, rus_prop, age, tenure, outsider, gov_opp, region_id, factor_year, out_prop) %>% na.omit()
fit6_multinom   <- multinom(gov_change_new3 ~ elite_consolidation + growth + grp_logged + unemp_delta + term_limits_dum + log(nat_res_2006 + 1) + log(distance_from_moscow_km + 1) +
              rus_prop + age + tenure + outsider + gov_opp + out_prop + factor_year, data = df_fit6)
df_fit7     <- suv_df %>% select(gov_change_new3, elite_consolidation,  growth, grp_logged, unemp_delta, term_limits_dum, ur_share, ur_share_2003, nat_res_2006, population_2018, distance_from_moscow_km, rus_prop, age, tenure, outsider, gov_opp, region_id, factor_year, out_prop) %>% na.omit()
fit7_multinom   <- multinom(gov_change_new3 ~ elite_consolidation + ur_share + ur_share_2003 + growth + grp_logged + unemp_delta + term_limits_dum +  + log(nat_res_2006 + 1) + log(distance_from_moscow_km + 1) +
              rus_prop + age + tenure + outsider + gov_opp + out_prop + factor_year, data = df_fit7)
df_fit8     <- suv_df %>% select(gov_change_new3, elite_consolidation,  growth, grp_logged, unemp_delta, term_limits_dum, ur_fed, ur_share_2003, nat_res_2006, population_2018, distance_from_moscow_km, rus_prop, age, tenure, outsider, gov_opp, region_id, factor_year, out_prop) %>% na.omit()
fit8_multinom   <- multinom(gov_change_new3 ~ elite_consolidation + ur_fed + ur_share_2003  + growth + grp_logged + unemp_delta + term_limits_dum +  + log(nat_res_2006 + 1) + log(distance_from_moscow_km + 1) +
              rus_prop + age + tenure + outsider + gov_opp + out_prop + factor_year, data = df_fit8)





se_fit1     <- se_return(fit1_multinom, df_fit1)
se_fit2     <- se_return(fit2_multinom, df_fit2)
se_fit3     <- se_return(fit3_multinom, df_fit3)
se_fit4     <- se_return(fit4_multinom, df_fit4)
se_fit5     <- se_return(fit5_multinom, df_fit5)
se_fit6     <- se_return(fit6_multinom, df_fit6)
se_fit7     <- se_return(fit7_multinom, df_fit7)
se_fit8     <- se_return(fit8_multinom, df_fit8)


stargazer(fit1_multinom, fit2_multinom, fit3_multinom, fit6_multinom, fit7_multinom, fit8_multinom, fit4_multinom, fit5_multinom,
          se = list(se_fit1, se_fit2, se_fit3, se_fit6, se_fit7, se_fit8, se_fit4, se_fit5), 
          omit = c('year'),
          add.lines = list(c('Observations', c(nrow(df_fit1), nrow(df_fit2), nrow(df_fit3), nrow(df_fit6), nrow(df_fit7), nrow(df_fit8), nrow(df_fit4), nrow(df_fit5))), 
                c('Appointments', c(nrow(df_fit1[df_fit1$gov_change_new3 == 'Local',]), nrow(df_fit1[df_fit1$gov_change_new3 == 'Outsider',]), nrow(df_fit2[df_fit2$gov_change_new3 == 'Local',]), nrow(df_fit2[df_fit2$gov_change_new3 == 'Outsider',]),
            nrow(df_fit3[df_fit3$gov_change_new3 == 'Local',]), nrow(df_fit3[df_fit3$gov_change_new3 == 'Outsider',]), nrow(df_fit6[df_fit6$gov_change_new3 == 'Local',]), nrow(df_fit6[df_fit6$gov_change_new3 == 'Outsider',]), 
            nrow(df_fit7[df_fit7$gov_change_new3 == 'Local',]), nrow(df_fit7[df_fit7$gov_change_new3 == 'Outsider',]), nrow(df_fit8[df_fit8$gov_change_new3 == 'Local',]), nrow(df_fit8[df_fit8$gov_change_new3 == 'Outsider',]),
            nrow(df_fit4[df_fit4$gov_change_new3 == 'Local',]), nrow(df_fit4[df_fit4$gov_change_new3 == 'Outsider',]), nrow(df_fit5[df_fit5$gov_change_new3 == 'Local',]), nrow(df_fit5[df_fit5$gov_change_new3 == 'Outsider',])
            ))), no.space = TRUE,
            covariate.labels = c('Civil Liberties', 'UR Vote Share in Regional Elections', 'UR Vote Share in State Duma Elections ', 'UR Vote Share in 2003 State Duma Election', 'Pluralism of Elites', 'Margins in Last Gubernatorial Elections (Before 2004)', 'Margins in Last Gubernatorial Elections (All)', 'Growth', 'GRP per capita (logged)', 'Unemployment', 'End of Term',
            'Natural Endowment (logged)', 'Distance from Moscow (logged)', 'Ethnic Russians (%)', 'Governor’s Age ', 'Tenure ', 'Incumbent Outsider Governor', 'Incumbent Opposition Governor', 'Prop. Outsider Gov.'))


pred_fit1       <- Effect('democracy_lower', fit1_multinom, vcov = mlogit.clust(fit1_multinom, df_fit1, 'region_id'), xlevels=list(democracy_lower = seq(5, 15, 1)))
pred_fit2       <- Effect('ur_share', fit2_multinom, vcov = mlogit.clust(fit2_multinom, df_fit2, 'region_id'), xlevels=list(ur_share = seq(10, 90, 2)))
pred_fit4       <- Effect('last_margin', fit4_multinom, vcov = mlogit.clust(fit4_multinom, df_fit4, 'region_id'), xlevels=list(last_margin = seq(0, 100, 2)))
pred_fit6       <- Effect('elite_consolidation', fit6_multinom, vcov = mlogit.clust(fit6_multinom, df_fit1, 'region_id'), xlevels=list(elite_consolidation = seq(4, 14, 1)))

Effect('ur_share', fit2_multinom, vcov = mlogit.clust(fit2_multinom, df_fit2, 'region_id'), xlevels=list(ur_share = quantile(df_fit2$ur_share, probs = c(.1, .5, .9))))
Effect('last_margin', fit4_multinom, vcov = mlogit.clust(fit4_multinom, df_fit4, 'region_id'), xlevels=list(last_margin = quantile(df_fit4$last_margin, probs = c(.1, .5, .9))))

pdf('pred_app_main_0521.pdf', width = 8, height = 8)
par(mfrow = c(2, 2))
plot(0, 0, xlim = c(5, 15), ylim = c(0, 25), xlab = 'Civil Liberty', ylab = 'Predicted Prob. of Appointments (%)', xaxt = 'n', yaxt = 'n', type = 'n', main = '(a) Model 1', bty = 'n')
lines(seq(5, 15, 1), pred_fit1$prob[, 'prob.Local'] * 100, lwd = 2, col = my_red)
lines(seq(5, 15, 1), pred_fit1$prob[, 'prob.Outsider'] * 100, lwd = 2, col = my_blue, lty = 2)
polygon(c(seq(5, 15, 1), rev(seq(5, 15, 1))), c(pred_fit1$upper.prob[, 'U.prob.Local'] * 100, rev(pred_fit1$lower.prob[, 'L.prob.Local'] * 100)),
            col = adjustcolor(my_red, alpha.f = .2), border = adjustcolor(my_red, alpha.f = .5))
polygon(c(seq(5, 15, 1), rev(seq(5, 15, 1))), c(pred_fit1$upper.prob[, 'U.prob.Outsider'] * 100, rev(pred_fit1$lower.prob[, 'L.prob.Outsider'] * 100)),
            col = adjustcolor(my_blue, alpha.f = .2), border = adjustcolor(my_blue, alpha.f = .5))
x_quantile      <- quantile(df_fit1$democracy_lower, probs = c(.1, .5, .9))
abline(v = x_quantile, lty = 3)
text(x = x_quantile, y = rep(18, 3), labels = c('10%', '50%', '90%'), pos = 4)
axis(side = 1, seq(5, 15, 2))
axis(side = 2, seq(0, 50, 5), las = 2)
legend('topleft', legend = c('Local Gov. App.', 'Outsider Gov. App.'), lty = 1:2, lwd = c(2, 2), col = c(my_red, my_blue), bty = 'n')


plot(0, 0, xlim = c(10, 90), ylim = c(0, 25), xlab = 'UR Share (%)', ylab = 'Predicted Prob. of Appointments (%)', xaxt = 'n', yaxt = 'n', type = 'n', main = '(b) Model 2', bty = 'n')
lines(seq(10, 90, 2), pred_fit2$prob[, 'prob.Local'] * 100, lwd = 2, col = my_red)
lines(seq(10, 90, 2), pred_fit2$prob[, 'prob.Outsider'] * 100, lwd = 2, col = my_blue, lty = 2)
polygon(c(seq(10, 90, 2), rev(seq(10, 90, 2))), c(pred_fit2$upper.prob[, 'U.prob.Local'] * 100, rev(pred_fit2$lower.prob[, 'L.prob.Local'] * 100)),
            col = adjustcolor(my_red, alpha.f = .2), border = adjustcolor(my_red, alpha.f = .5))
polygon(c(seq(10, 90, 2), rev(seq(10, 90, 2))), c(pred_fit2$upper.prob[, 'U.prob.Outsider'] * 100, rev(pred_fit2$lower.prob[, 'L.prob.Outsider'] * 100)),
            col = adjustcolor(my_blue, alpha.f = .2), border = adjustcolor(my_blue, alpha.f = .5))
x_quantile      <- quantile(df_fit2$ur_share, probs = c(.1, .5, .9))
abline(v = x_quantile, lty = 3)
text(x = x_quantile, y = rep(18, 3), labels = c('10%', '50%', '90%'), pos = 4)
axis(side = 1, seq(10, 100, 10))
axis(side = 2, seq(0, 50, 5), las = 2)
legend('topleft', legend = c('Local Gov. App.', 'Outsider Gov. App.'), lty = 1:2, lwd = c(2, 2), col = c(my_red, my_blue), bty = 'n')
temp    <- suv_df %>% filter(gov_change == 1 & is.na(ur_share) == FALSE)
for(i in 1:nrow(temp)){
        temp_col        <- ifelse(temp$out_change[i] == 1, my_blue, my_red)
        temp_lty        <- ifelse(temp$out_change[i] == 1, 2, 1)
        lines(rep(temp$ur_share[i], 2), c(0, -5), lty = temp_lty, col = temp_col)
}

plot(0, 0, xlim = c(4, 14), ylim = c(0, 25), xlab = 'Elite Consolidation', ylab = 'Predicted Prob. of Appointments (%)', xaxt = 'n', yaxt = 'n', type = 'n', main = '(d) Model 4', bty = 'n')
lines(seq(4, 14, 1), pred_fit6$prob[, 'prob.Local'] * 100, lwd = 2, col = my_red)
lines(seq(4, 14, 1), pred_fit6$prob[, 'prob.Outsider'] * 100, lwd = 2, col = my_blue, lty = 2)
polygon(c(seq(4, 14, 1), rev(seq(4, 14, 1))), c(pred_fit6$upper.prob[, 'U.prob.Local'] * 100, rev(pred_fit6$lower.prob[, 'L.prob.Local'] * 100)),
            col = adjustcolor(my_red, alpha.f = .2), border = adjustcolor(my_red, alpha.f = .5))
polygon(c(seq(4, 14, 1), rev(seq(4, 14, 1))), c(pred_fit6$upper.prob[, 'U.prob.Outsider'] * 100, rev(pred_fit6$lower.prob[, 'L.prob.Outsider'] * 100)),
            col = adjustcolor(my_blue, alpha.f = .2), border = adjustcolor(my_blue, alpha.f = .5))
x_quantile      <- quantile(df_fit6$elite_consolidation, probs = c(.1, .5, .9))
abline(v = x_quantile, lty = 3)
text(x = x_quantile, y = rep(18, 3), labels = c('10%', '50%', '90%'), pos = 4)
axis(side = 1, seq(4, 14, 2))
axis(side = 2, seq(0, 50, 5), las = 2)
legend('topleft', legend = c('Local Gov. App.', 'Outsider Gov. App.'), lty = 1:2, lwd = c(2, 2), col = c(my_red, my_blue), bty = 'n')


plot(0, 0, xlim = c(0, 100), ylim = c(0, 25), xlab = 'Margins in Last Gubernatorial Elections (%)', ylab = 'Predicted Prob. of Appointments (%)', xaxt = 'n', yaxt = 'n', type = 'n', main = '(c) Model 7', bty = 'n')
lines(seq(0, 100, 2) , pred_fit4$prob[, 'prob.Local'] * 100, lwd = 2, col = my_red)
lines(seq(0, 100, 2) , pred_fit4$prob[, 'prob.Outsider'] * 100, lwd = 2, col = my_blue, lty = 2)
polygon(c(seq(0, 100, 2) , rev(seq(0, 100, 2) )), c(pred_fit4$upper.prob[, 'U.prob.Local'] * 100, rev(pred_fit4$lower.prob[, 'L.prob.Local'] * 100)),
            col = adjustcolor(my_red, alpha.f = .2), border = adjustcolor(my_red, alpha.f = .5))
polygon(c(seq(0, 100, 2) , rev(seq(0, 100, 2) )), c(pred_fit4$upper.prob[, 'U.prob.Outsider'] * 100, rev(pred_fit4$lower.prob[, 'L.prob.Outsider'] * 100)),
            col = adjustcolor(my_blue, alpha.f = .2), border = adjustcolor(my_blue, alpha.f = .5))
x_quantile      <- quantile(df_fit4$last_margin, probs = c(.1, .5, .9)) 
abline(v = x_quantile, lty = 3)
text(x = x_quantile, y = rep(18, 3), labels = c('10%', '50%', '90%'), pos = 4)
axis(side = 1, seq(0, 100, 10))
axis(side = 2, seq(0, 50, 5), las = 2)
legend('topleft', legend = c('Local Gov. App.', 'Outsider Gov. App.'), lty = 1:2, lwd = c(2, 2), col = c(my_red, my_blue), bty = 'n')
temp    <- suv_df %>% filter(gov_change == 1 & is.na(last_margin) == FALSE)
for(i in 1:nrow(temp)){
        temp_col        <- ifelse(temp$out_change[i] == 1, my_blue, my_red)
        temp_lty        <- ifelse(temp$out_change[i] == 1, 2, 1)
        lines(rep(temp$last_margin[i] , 2), c(0, -5), lty = temp_lty, col = temp_col)
}
dev.off()



